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ABSTRACT 

Context. There are numerous examples of accretion discs in binary systems where the disc midplane is believed to be inclined relative to the 
binary orbit plane. 

Aims. We aim to examine the detailed disc structure that arises in a misaligned binary system as a function of the disc aspect ratio h, viscosity 
parameter a, disc outer radius R, and binary inclination angle y F . We also aim to examine the conditions that lead to an inclined disc being 
disrupted by strong differential precession. 

Methods. We use a grid-based hydrodynamic code to perform 3D simulations. This code has a relatively low numerical viscosity compared 
with the SPH schemes that have been used previously to study inclined discs. This allows the influence of viscosity on the disc evolution to be 
tightly controlled. 

Results. We find that for thick discs (h = 0.05) with low a, efficient warp communication in the discs allows them to precess as rigid bodies 
with very little warping or twisting. Such discs are observed to align with the binary orbit plane on the viscous evolution time. Thinner discs 
with higher viscosity, in which warp communication is less efficient, develop significant twists before achieving a state of rigid-body precession. 
Under the most extreme conditions we consider (h = 0.01, a = 5 x 10 -3 and a = 0.1), we find that discs can become broken or disrupted by 
strong differential precession. Discs that become highly twisted are observed to align with the binary orbit plane on timescales much shorter 
than the viscous timescale, possibly on the precession time. 

Conclusions. We find agreement with previous studies that show that thick discs with low viscosity experience mild warping and precess 
rigidly. We also find that as h is decreased substantially, discs may be disrupted by strong differential precession, but for disc thicknesses that 
are significantly less (h = 0.01) than those found in previous studies (h = 0.03). 



1. Introduction 

Accretion discs are found in a number of astrophysical environ- 
ments. These include protostellar discs around T Tauri stars, 
and around compact objects in cataclysmic variable and X- 
ray binary systems. Most T Tauri stars associated with low- 
mass star-forming regions are observed to be members of bi- 
nary systems, where the peak in the d istribution of separa - 



tions occurs at approximately 30 AU dLeinert et al. 1 11993). 



Studies of the tidal interaction between a circumstellar disc 
and a binary companion indicate that tidal truncation oc- 
curs at a ratio of binary sepa ration to disc radius D/R - 



3 for an equal ma s s binary ( Artvmo wicz & Lubow 1 1 1994 



Larwood et all 19961: iKley & Nelsonll2008h . Given that 30 AU 
is smaller than typica l disc sizes observed in T Tauri systems 



(Ed wards et al. II 19871) . protostellar discs in many binary sys- 



tems are often going to experience strong tidal effects. 

There are observational indications that the disc and binary 
orbital plane are not always coplanar. The binary system HK 
Tau has been imaged, and shows compelling evidence for a 
disc whose midplane is misaligned with the orbit of the bi- 



nary compan ion (IStapelfeldt et alJll998h . iTerquem & Bertout 
dl993lll996h ', tave suggested that the reprocessing of radiation 
from the central star by a tilted and warped accretion disc could 
account for the high spectral index of some T Tauri stars. A 
number of precessing jets have been observed in star forming 
regions, and the jet precession may indicate a n underlying pre 



cessing disc from which th e jet is launched dKonigl & Ruden 



1993 : Terquem et all 19991) . In the case of close X-ray binaries, 
the objects SS433 and Her X- 1 show evidence of precessing 
jets and/or accretion discs dBovnton et al. Il98fj lKatz| ["l980: 
Margonll984l:IPetterson Ill977l) . 



An early study of the equations describing the be haviour of 
a thin non-coplanar viscous disc was undertaken by lPettersonl 
dl977l) . The resulting equations, however, lacked the property 
of conserving global an gular momentum. This was a ddressed 
in subsequent work by IPapaloizou & Pringld ( 11983b . where 



they showed that warps in discs for which h < a < 1 evolve 
diffusively on a timescale shorter than the usual viscous ac- 
cretion by a factor or 2 , where a is the dimensionless viscosity 
coefficient ( Shakura & Sunvaevlll973 ). and h is the disc aspect 
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ratio H/r. The oppos i te reg ime a < h has been investigated 
by iPapaloizou & Linl (119951) using linear perturbation theory. 
They demonstrated that the governing equation for the disc tilt 
changes from being a diffusion equation to a wave equation in 
this physical regime, showing that warps propagate as bending 
waves with a speed corresponding to half the sound spe ed, c s /2. 
An an alytic study of warped discs was undertaken by lOgilvie 
(2000) in which the mildly non linear evolution of disc warps 
was examined. 

The tidal perturbation of an inviscid disc by a com- 
panion on an inclined circ ular orbit has been studied by 
Papa loizou & T erquem d 1995b using linear theory. They inves- 
tigate the disc response to tidal perturbations that have odd 
symmetry with respect to the disc midplane and consider both 
zero and non-zero perturbing frequencies. They show that the 
zero-frequency (secular) mode leads to rigid precession of the 
disc if the sound crossing time through the disc is smaller than 
the differential precession time. Furthermore, they estimate that 
the timescale for damping the inclination is of the same order as 
the accretion timescale of the disc. A study of the twisting and 
warping of viscous discs in binary systems was undertaken by 
Lubow & Ogilvid pOOO) using linear theory, and it was shown 
that discs whose outer disc radii are smaller than the truncation 
radius maybe unstable to tilting. 

Numerical simulations of tidally interacting discs which 
are not coplanar with the binary orbit have been performed by 
Larwood et al.1 (119961) using SPH simulations. They found that 
the disc precesses approximately as a rigid body as long as the 
disc aspect ratio is large enough. For smaller aspect ratios their 
results suggest that the disc develops a modest warp, while for 
h < 0.03 the disc becomes disrupted as a consequence of dif- 
ferential precession. Their results suggest that a disc can split- 
up into two distinct parts that behave independently. These 
SPH simulations also showed that the inclination evolves on 
the vi scous timescale, as expected from lPapaloizou & Terqueml 
dl995l) . 

The goal of the present study is to examine in detail the 
structure of misaligned accretion discs in close binary systems 
as a function of the important physical parameters: h, a, the 
outer disc radius R, and the inclination angle jf. The binary 
companion and the central star are assumed to be of equal mass. 
We use a grid-based code (NIRVANA) to perform this study, 
and this code has a relatively small numerical viscosity com- 
pared with the SPH schemes used to perform previous studies. 
This allows us to tightly control the disc viscosity, and therefore 
examine in more detail its influence on the disc evolution. We 
consider a broad range of disc parameters in which warps prop- 
agate either as bending waves h > a or diffusively h < a. As 
well as examining the quasi-steady disc structures which arise 
because of tidal interaction with the inclined companion, we 
also examine the conditions under which a disc becomes dis- 
rupted due to strong differential precession. In particular we are 
interested in examining whether the following criterion holds: a 
disc will achieve a state of rigid-body precession if warp prop- 
agation across the disc (either through bending waves or diffu- 
sion) occurs on a timescale shorter than the differential preces- 
sion time. 



The plan of the paper is as follows. In Sect. 2 we present 
the basic equations of the problem. In Sect. 3 we describe the 
numerical methods used in the simulations, and in Sect. 4 we 
discuss the linear theory of bending waves and calibrate the 
simulation code against calculations based on linear theory. In 
Sect. 5 we present our results for the tilted, tidally interacting 
discs, and in Sect. 6 we discuss our results and present our con- 
clusions. 



2. Basic equations 

We consider the evolution of a thin, gaseous, viscous disc in 
orbit around a central star, where the disc is perturbed by a 
companion star orbiting in a plane which is not coincident with 
the disc midplane. For a range of physical parameters, it is ex- 
pected that the disc will precess rigidly around the angular mo- 
mentum vector of the binary system. The equations of conti- 
nuity and motion for a viscous fluid in a precessing reference 
frame may be written as: 



— + pV.v = 
Dt H 



(1) 



— = —VP - + S w «. - 2Q xv-fix(Jlxr) (2) 
Dt p 

where 
D d 

— = — + v V 

Dt dt 

is the convective derivative, p is the density, v the velocity and P 
the pressure. The precession frequency is given by |£2|, and the 
disc angular momentum vector is assumed to precess around a 
vector pointing in the direction of £1. The gravitational poten- 
tial is We adopt a locally isothermal equation of state: 



P = c s (rfp 



(3) 



where the isothermal sound speed is defined by c s (r) = h ■ VKep- 
Here h is the aspect ratio of the disc, H/r, and Vjtep is the local 
keplerian velocity. S V!JC is the viscous force per unit mass: 



Snsc = -vt 
p 

v, Ik' iv T is the viscosity stress tensor: 
tdvi dvj 2 \ 

r -^ + ^-3*H 



(4) 



(5) 



We adopt the standard 'alpha' model (Shakura & Sunyaev 
1973) to specify the kinematic viscosity v = ac s H. We solve 
the above equations numerically, using a system of spherical 
coordinates r = (r, 6, <p). 

2.1. Orbital configuration 

We work in a reference frame in which the origin of the coordi- 
nate system is fixed on the central star, and the secondary star 
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moves on a circular orbit about the origin with position vec- 
tor D. The gravitational potential, *P, at any position vector r is 
therefore given by: 



GM P GM S GM S r ■ D 



D| 



D 3 



(6) 



where G is the gravitational constant, and Mp and Ms are the 
primary and secondary masses, respectively. The first term is 
due to the primary star, and the remaining two terms are due 
to the companion star. The last indirect term accounts for the 
acceleration of the origin of the coordinate system, which co- 
incides with the location of the primary star. Note that there are 
no contributions to the potential due to the disc since we do not 
include disc self-gravity in our calculations, and the stars do not 
experience a gravitational force due to the disc. Consequently 
the binary orbital parameters remain unchanged. 

Hence the secondary star feels the acceleration due to the 
central star and the indirect component of the potential. In a 
non-precessing reference frame, centred on the central star, 
whose x — y plane coincides with the orbital plane of the bi- 
nary, the position vector of the companion in a circular orbit 
may be written as: 

D = D ■ (cos (wf)ei + sin (a>t)t2), (7) 

where D is the constant binary separation and u> = 
^G(Mp + Ms)/D 3 is the binary orbital angular velocity. 



3. Numerical methods 

The system of equations described in Sect. [2] is inte- 
grated using the grid-b ased hydrodynamics code NIRVANA 
( Zie gler&Yorkel 1997b . adapted to solve the equations in a pre- 
cessing reference frame. This code uses operator-splitting, and 
the advection routin e uses a second -order accurate monotonic 
transport algorithm ( van Leerlll977 ). 



3.1. Boundary conditions 

Periodic boundary conditions were applied in the azimuthal 
direction. At all other boundaries outflow conditions were 
adopted using zero-gradient extrapolation. Velocities at the ra- 
dial boundaries were set to ensure that zero viscous stress oc- 
curs there. 



3.2. Units 

We adopt a system of units such that the unit of mass is equal 
to that of the central star, Mp, the unit of radius is equal to 
the radius of the inner boundary of our computational domain, 
and we set the gravitational constant G — 1 . The unit of time 
then becomes 1 /Q/K 1), where Q# is the keplerian angular ve- 
locity evaluated at r = I. When discussing simulation results, 
however, we will express time in units of the orbital period 
at r — 10, the nominal outer disc radius in the simulations. 
Thus we describe a time interval of 27iQ. K {\0y l as 'one orbit'. 
Inclination and precession angles are described in units of de- 
grees. 



3.3. Reference frames and initial conditions 

The numerical domain extends radially from r = [1,12], and 
azimuthally from (p = [0, 360°]. The outer radius, R, may differ 
from simulation to simulation. The meridional domain covers 
the range = [90° - AO, 90° + AO], where again AO may vary 
between the simulations. 

We perform simulations in two different reference frames. 
For models in which the disc is expected to develop rigid pre- 
cession, we solve the equations in a frame which precesses 
around the angular momentum vector of the binary. Given a 
sensible choice of the precession rate, this approach ensures 
that the disc midplane always stays close to the equatorial plane 
of the computational domain where = 90°, and allows sim- 
ulations to be performed with large inclinations of the binary 
orbit. If evolved in a non-precessing frame, the disc would pre- 
cess around the binary angular momentum vector, causing it 
to eventually intersect with the upper and lower meridional 
boundaries of the computational domain if the binary inclina- 
tion angle is large. Adopting a precessing reference frame al- 
lows simulations to be conducted with relatively small merid- 
ional domains, even for large binary inclinations, thus substan- 
tially reducing the computational expense involved. 

For models in which substantial differential precession of 
the disc is expected to arise, adopting a frame with a single pre- 
cession frequency does not solve the problem described above, 
since some parts of the disc inevitably intersect the meridional 
boundaries. As these models are of significant interest from 
the astrophysical point of view (they will tend to be thin discs 
with large viscosity, similar to those which occur in X-ray bina- 
ries for example), we have computed them in a non precessing 
frame, but have been forced to use large meridional domains 
and smaller binary inclination angles. 

We denote coordinates in the precessing frame as f, while 
coordinates in the non-precessing binary frame are denoted r. 
The coordinates in the code frame are denoted with rc, where 
the code frame can be one of either the precessing frame or 
the binary frame. The transformation of any vector e from the 
binary into the precessing frame is given by: 

e = R z (Cl F t)R x (j F )e (8) 

with rotation matrices: 

'10 

Rx(Jf) = cos(y f ) sin(y F ) 

k - sin(y F ) cos(y F ) J 

'cos(Q F f) -sin(Q F f) 0' 
R z (Q. F t) = sin(Q F f) cos(fi F f) . (9) 
{ lj 

Thus to transform a vector from the binary into the precessing 
frame, we first rotate around the x-axis by an inclination an- 
gle yp and then rotate the resulting vector around the z-axis by 
an precession angle Q. F t. Here Q. F is the precession rate of the 
precessing frame and t is the time elapsed. This transformation 
is particularly useful when setting up the initial conditions for 
the simulations which are performed in the non-precessing bi- 
nary frame, since we assume that the binary orbit lies in a plane 
which is coincident with the equatorial plane of the computa- 
tional grid, with the disc being inclined by an angle yp. 
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3.3.1 . Model set-up in precessing frame 



3.3.2. Model set-up in binary frame 



In the precessing frame the code midplane (defined to reside 
where 9c — 90°) coincides with the disc midplane initially, so 
that (r c = f). The initial density profile is chosen to satisfy hy- 
drostatic equilibrium in the ^-direction, which yields the result: 



p(r, 9) = r ( [sin (9)] 



(£-<f +1 >) 



(10) 



where ( = | is the radial power law exponent of the density 
profile. Note that this expression has the limit: 



lim p(r, 9) — r ( exp 
e-»o 



2h 2 



(11) 



which is the familiar form of the density profile for thin discs. 
The velocity in the direction is chosen to balance the forces 
in the radial direction: 



v 2 = l -{\-^+\)h 2 ). 



(12) 



Consequently, the azimuthal velocity profile is slightly subke- 
plerian because of the radial pressure gradient. The velocities 
in the r and 9 direction are chosen to be zero initially. 
We set £1 in Eq. Q to point in the same direction as the binary 
angular momentum vector: 



= Qf-e3 = Q.f(- sin(yf)e2 + cos (7^)63) 



where the last equality is derived using the transformation 
given by Eq. (0. The precession rate may be estimated using 
the expression: 



-I 



3GMs 



cos (y F ), 



(14) 



where X = f p r d9 is the disc surface density. Hence the pre- 
cession frequency is determined by the ratio of the integrated 
external torque exerted by the second ary to the integrated an- 
gular momentum content of the disc (IPapaloizou & Terquem 
Il995l) . For the velocity and density profiles given by Eqs. (fT2b 
and ( TTOb . respectively, this evaluates to: 



3 R 3 M s 



Q d (R) 1 D 3 M P 



cos (y F ), 



where Q.j(R) is the orbital angular velocity of the disc at radius 
R. The binary orbit given by Eq. (0, expressed in precessing 
frame coordinates, is given by: 



D = D 



cos ([to - Qf]f) 
sin ([to - Qf]f) cos(y F ) 
sin ([to - Of]?)) sin (Jf) , 



(16) 



Thus an observer moving with the disc sees an increased binary 
frequency to - Q.p due to the retrograde precession of the disc 
(i.e. Q F is negative). 



When working in the non-precessing binary frame we set the 
precession frequency £1 = in Eq. ©. Since the equatorial 
plane of the computational domain coincides with the binary 
orbit plane (rc = r), the disc is inclined with respect to the 
equatorial plane of the computational grid by an inclination an- 
gle jf. We can use the transformation defined by Eq. © at 
t — to relate the polar angles of the precessing frame to the 
polar angles of the binary frame. One particularly useful rela- 
tion is given by: 



sin 2 (9) = cos 2 {(/)) sin 2 (9) + cos 2 (y F ) sin 2 (0) sin 2 (9) 
+ sin 2 (y F ) cos 2 (9) 

- 2 cos 2 (yp) sin (yp) sin (9) cos (9) sin (cf>) 



(17) 



which can be used to determine the initial density profile given 
by Eq. ( TTOb . The initial velocity is given by: 

= v (-sin(0)ei +cos(0)e 2 ) 

= v^(- sin (0)d + cos ($) cos (jf)^2 - cos 0) sin (7^)63) 

Using the relationship between the polar angles one can find 
an expression for that is written entirely in terms of binary 
frame angles. The different components in the radial, azimuthal 
and meridional direction are then given by: 







(13) = 



cos (yp) sin (9) - sin (<p) sin (yp ) cos (9) 



sin (9) 



sin (yp) cos (tp) 
sin (9) 



(18) 



where v$ is given by Eq. (fT2l and sin (9) by Eq. ( TT7l >. It can 
be shown that this set of profiles satisfies the Euler equations 
given by Eqs. ([T]) and ® with S,, L , C = 0, M s = 0, = and 
yp - constant, since in a spherically symmetric potential there 
is no preferred direction for the disc rotation axis. 

3.4. Definition of precession and inclination angles 

In order to follow the evolution of the disc structure we calcu- 
late the total angular momentum vector, L(r), for disc material 
contained in each spherical shell in the numerical domain: 



(15) L(r) 



pin pj+AS 
I d<P \ 1 

Jo Jf-Ae 



I r 2 sin(9)d96r 



(19) 



where 6r is the radial grid spacing and 1 is the angular momen- 
tum density: 



1 = pr x v = p(-ry^6 c + rv g (f> c ). 



(20) 



9c and 4>c are tne code unit vectors in the meridional and az- 
imuthal directions, respectively. Note that when working in the 
precessing frame a term involving the precessional velocities 
should be added to this expression. In practice we find that this 
is of negligible magnitude, because of the slow rate of preces- 
sion, and therefore the term has been omitted. For disc material 



M.M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 



5 



located in a given spherical shell we locate the inclination an- 
gle, 5, (equal to y F at time t = 0) between the disc and binary 
angular momentum vectors through: 



cos (6) — 



UbIILI 



(21) 



where Jb represents the binary angular momentum vector. The 
precession angle, /3, measured in the binary plane can be de- 
fined through: 



cos(/3) 



(Jb x L).u 
|J B xL| ' 



(22) 



where u is an arbitrary reference unit vector in the binary plane. 
We choose u = — ei , so that the initial precession angle (3 — by 
this definition. When working in the precessing frame -ei has 
to be expressed in precessing frame coordinates and becomes 
time dependent. 

When discussing the results of the simulations we use the 
the following terminology: a warp occurs if the angle 5 be- 
comes a function of radius in the disc; a disc develops a twist 
when f3 becomes a function of radius; a disc is said to be bro- 
ken if either j3 or 5 change discontinuously at some radius such 
that the disc separates into two independently precessing parts; 
a disc is said to be disrupted if f3 varies smoothly by more than 
180 degrees across the disc radius because of differential pre- 
cession. 

3.5. Calculation of torque contributions 

Later in this paper we will be interested in measuring the dif- 
ferent contributions to the changes in the precession and incli- 
nation angles for disc material located at different radii in the 
computational domain. As a first step we calculate the rate of 
change of the angular momentum vector L(r) associated with 
disc material located at radius r, which involves integrating 
over the individual spherical shells in the computational do- 
main. 

Using the continuity and momentum Eqs. (Q~|) and (f2]i to 
express the velocity and density changes in Eq. ( 1201 ). and as- 
suming that the density vanishes at the meridional boundaries, 
we find the result: 

L(r)4 f W + LV) + L s W (23) 
where a dot denotes a time derivative. We have that 

dip ' d0sm{6)\r 2 -\-vA J (24) 

is the change due to the divergence of a radial angular momen- 
tum flux. The change due to viscous interaction with neighbor- 
ing shells is given by: 

i/(r)= d<f> r 3 sin (0) d0\-T r(P c + T rg <p c )\ ;(25) 

Jo Jf-A6> a 

Thus the z-component of the angular momentum vector can 
only be changed by the T r $ component of the viscous stress ten- 
sor, while the x and y components can also be changed due to 
radial shear of vertical motion via T r g. The last term in Eq. ( 1231 



which induces a change of the angular momentum vector of 
disc material contained in a spherical shell is due to the gravi- 
tational interaction with the secondary star: 



lV) = lV) + lV), 



(26) 



where the angular momentum change due to torques is given 
by: 



>0 Jf-Afl 

and the torque density t is given by: 



nix p%+A9 

L r (r) = I d(f> I ' t-r 2 sin(0)d06r 

Jo J?-A9 



-pr x V*F = p 



1 ff¥ 
sin (8) Ihj) 



~do ( ' 



(27) 



(28) 



When working in the precessing frame there is an additional 
term causing angular momentum change due to Coriolis and 
centrifugal forces: 



t A (r) = I dip \ pr 2 sin (0) d6 (-rAv^c + rAv e tp c ) 8r 
Jo J | -AS 

with 



Av = (-2Q x v - Q x (a x r)) <p c 
Av g = (-2Q x v - x (« x r)) C . 



(29) 



In the following we are only interested in the torque compo- 
nents expressed in precessing frame coordinates. Thus all the 
torque components have been calculated now for simulations 
performed in the precessing frame. However, when working in 
the binary frame, we perform the above calculations and then 
project each of the torque components onto precessing frame 
coordinate axes to ensure a uniform approach to monitoring 
the torque contributions. Since the latter are time dependent, 
the additional term when working in the binary frame is given 
by: 

t A (r) = -£l X L 

This term accounts for the angular momentum change mea- 
sured in the precessing frame that arise because of the preces- 
sion of the frame. 



3.6. Calculation of precession and inclination rates 

We can now relate the above torques to the rate of change of the 
precession and inclination angles. When working in the pre- 
cessing frame, the inclination angle 6 and precession angle /3 
for disc material in a given radial shell, as defined by Eqs. < f2Tb 
and d22l . are given by: 



cos(j6) 



cos(5) = 



cos(£2p f)(sin(y/r)Lz + cos(yf )Ly) + sin(Q.ft)Lx 

[(sin(y F )L z + cos(y f )L Y ) 2 + L 2 X ] 2 
■ sin(y F )L F + cos(y F )L z 



|L| 



(30) 



Note that the vector components L x , L Y and L z in the above 
expression, and in the remaining expressions in this subsec- 
tion, should have hat symbols because they are measured in 
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the precessing frame. We neglect these, however, in order to 
simplify the notation. For sufficiently rigidly precessing discs 
the deviation from the mean precession frequency Q. F will be 
small, and the disc midplane will remain close to the equatorial 
plane of the computational domain. In this case we can use the 
approximation: 



(31) 



With this approximation, Eq. ( f30b gives to first order: 
Ly 



JF + 



1 L x 
sin (y F ) L z ' 



(32) 



Thus for sufficiently rigidly precessing discs the total inclina- 
tion is the sum of the inclination of the precessing frame and a 
small deviation term which is proportional to the 3'-component 
of the angular momentum vector, measured in the precessing 
frame. Likewise the precession angle is the sum of precession 
due to the frame and a small deviation term that is proportional 
to the jc-component of the angular momentum vector. The ad- 
vantage of this approximation is that the precession and inclina- 
tion angles become linear functions of the torque components. 
This allows us to write: 



assuming an isotropic viscosity. iPapaloizou & Lira ( 11995b in- 
vestigated the regime in which a < h, and showed under that, 
for low frequency modes, the governing equation for the disc 
tilt changes from being a diffusion equation to a wave equation, 
such that warps are communicated through the disc via propa- 
gation of bending waves. In a keplerian disc these are expected 
to pr opagate with little dispersio n at half the sound speed, c s /2. 

INelson & Papaloizoul (|l999h used smooth particle hydro- 
dynamic simulations to examine the propagation of free bend- 
ing waves, comparing their numerical results with linear cal- 
culations. Overall, they found that the SPH code managed to 
model the propagation of bending waves quite accurately as 
long as the warp amplitude remains small such that linear the- 
ory holds. One drawback associated with using SPH for that 
problem, however, was the fairly large numerical viscosity as- 
sociated with the scheme which is difficult to control and spec- 
ify ab initio. 

In this section we compare our numerical results with 
a linear calculation, adopting a sim ilar set-up to the one 
used by INelson & Papaloizo ul (ll999h . The equations describ- 
in g the evolution of linear ben ding waves have been derived 
by Nelson & P apaloizou ( 1999h . We do not reproduce the full 
derivation here, but rather mention some of the salient points 
associated with the derivation. 



86 _ I86Y ld6Y Id6 
dt ~ [dtj + \8tj + \8t 

dt \ dt) [dt \dt 



(33) 



4.1. An initial value problem 

Small-amplitude warps can be considered to be linear pertur- 
bations of a disc with midplane initially coincident with the 
(ei, ex) plane. Taking the <p dependence of the perturbations 
where each of the rates are calculated according to Eq. d32]): to be cc e'"^, with m = 1 for global warps, the Euler equa- 
tions written in cylindrical coordinates (r, <p, z) can be linearised 
dPapaloizou & Linll 19951: INelson & Papaloizoull 1999b . An ana- 
lytical solution to the linearised equations corresponds to the 
case of a rigid tilt. In this case the velocity and pressure pertur- 
bations are given by: 



(i) 

ld£\ K 
\dt 



- 1 j'K _ Ly TK 

" L Z Y Ll z 



1 



sin(y F ) 



Lz Lx 



hi j'k 

Ll z 



(34) 



Hence, in circumstances where the approximation expressed 
in Eq. OTb holds, we can measure separate contributions to 
the inclination and precession rate coming from the radial flux 
(Eq. l24l . viscous stress (Eq. EBT l and gravitational interaction 
with the secondary star (Eq. |26*| |. When discussing contribu- 
tions to differential precession rates later in the paper, we will 
omit the constant mean precession rate D. F of the precessing 
frame in the expression for the total rate given by Eq. ( T33] >. 

4. Linear theory of free warps and bending waves 

In order to calibrate the code, we now examine how well 
it is able to model the propagation of free bending waves, 
and compare the results with linear theory. The linear the- 
ory of_war£S_jn_accretiondiscs was initially investigated 
by IPapaloizou & Pringld dl983b . Their analysis is valid in 
a regime in which the dim ensionless viscosity parameter a 
dShakura & Sunyaevlll973l) lies in the range h < a < 1. In 
this regime, globally warped discs were found to evolve dif- 
fusively, with a diffusion coefficient larger than that associ- 
ated with mass flow through the disc by a factor of ~ ^j, 



v' r =-zCl K d 



v '<p =-izS- 
v' z = r£l K 8 

cl 

W =p' — 
P 



dr 



-irzQ K 8 



(35) 



where Q.k is the Keplerian angular frequency and 5 is a small 
constant inclination. For large scale warps, it is expected that 
the local inclination 6 varies on a length scale significantly 
larger than the disc thickness H = h ■ r. It is then reasonable 
to assume the inclination is also approximately independent of 
z. Using these assumptions, and the rigid tilt solution d35b . the 
z dependence of the linearised equations disappears when in- 
troducing the new quantiti es q' r _ ^ and q'^ _ 4- . The re sulting 
expressions are given by INelson & Papalo izou (1 19991) in the 



form: 



— + iLl K q r - 2il K 



— i- 



d(rCl 2 K 6) 
dr 



dt 



, , 1 d(r 2 Cl K ) 2 
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Fig. 1. Time evolution of free warp for h=0.03, 6 



MAX 



5°. NIRVANA results (blue solid line); linear results (red dashed line). 



dv' , , 

— + in K v, = irQld 
dt z K 

JFr Q* — + i£l K 5 = + + ;2v, 

V dt / r or r * 



(36) 



where 



T 



«_/ —ex 
«7 —a 



^2 



pz 2 dz 



Equations d36i > are a set of coupled differential equations which 
describe the evolution of the inclination 6 as a function of r and 
t. Thus they can be used to follow the evolution of a bending 
disturbance from initial data. Below we use the time-dependent 
solution calculated in this way to test NIRVANA'S ability to 
propagate bending disturbances. 

Assuming a time dependence of the perturbations oc e'°~', 
with the wave frequency obeying the relation <x <sc Q. K for 
slowly varying warps, one can derive algebraic equations for 
the perturbed quantities. These show that warps induce hor- 
izontal motions and vertical shear in the disc through terms 

proportional to , ^ j ( Papaloizou & Linlll995 ). Physically, 

these horizontal motions are generated by pressure forces in the 
disc which arise from the vertical misalignment of neighboring 
regions of the disc midplane. They are responsible for driving 
the bending wave, which propagates with approximately half 
the sound speed. 

The main effect of viscosity in the disc is to damp these hor- 
izontal motions, reducing the efficacy of bending wave prop- 



agation. As the viscosity increases and a > h, the bend- 
ing disturbance begins to evolve diffusively on a time scal e 
that scales inversely with v/(2a 2 ) JPapaloizou & Pringlj 1983b . 
where this latter quantity is now the diffusion coefficient asso- 
ciated with warp propagation in this regime. If the viscosity is 
small enough that its effect on the unperturbed quantities can 
be neglected, then the set of equations d36l > can be extended 
to include the effect of a small viscosity by the replacement 
dPapaloizou & Linlll995l) : 



d 



, d 

-> I — + (z + a)Q K 



(37) 



If the viscosity is large, however, then the complete viscous 
stress tensor should be included in the equations of motion, 
making the analysis rather complicated, since there would then 
be accretion velocities in the unperturbed state. 

4.2. Setting up initial data 

A number of simulations using NIRVANA have been per- 
formed, in which warps with different amplitudes in discs with 
different parameters were set up and allowed to evolve. We 
compare these solutions with those obtained using the linear 
theory expressed in Eqs. (f36t using the same set of initial con- 
ditions. The linear calculations employed a one dimensional fi- 
nite difference scheme for which we used 1000 equally spaced 
grid cells. Starting with a gaussian initial inclination profile: 



6(r) = 6 M Axe 



(-k-r M/0 | 2 ) 
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h=0.01, <5„«=1.67 h = 0.03, <5„„ = 2 




Fig. 2. Final stage of evolution for different sets of disc parameters. NIRVANA solution is shown by the solid blue line, the linear 
calculation by the red dashed line. 



with rMiD — 5 to be centred in the middle of a disc with in- 
ner radius r — 1 and outer radius r — 9, the set of equations 
(f36t was integrated forward in time. We chose the viscosity 
in these test simulations to be a = 10~ 3 <K h since we are 
mainly interested in the regime of propagating bending waves, 
and because the linear solution presented here becomes invalid 
for large viscosities. The NIRVANA code was set up in an in- 
ertial frame using the expressions ( fTUb and ( [T8l to specify the 
initial density and velocity, in combination with Eq. ( fTTI ). with 
7fW = 8(f) now being the disc tilt which is a function of r. In 
these simulations we used N r = 404 cell in radius, 5 cells per 
pressure scale height in the meridional direction (and different 
numbers of vertical scale heights, depending on the initial warp 
amplitude, as described below) and N$ = 300 cells in azimuth. 
We used the same boundary condition as for the global simula- 
tions, which were described in section 3.1. 

4.3. Results 

An example result is shown in FigQ] which has a disc as- 
pect ratio of h — 0.03 and maximum initial inclination angle 
8 max - 5°. In this particular simulation the meridional domain 
contains 15 pressure scale heights (A6 = 12.89°) in total. FigQ] 
shows the inclination as a function of radius at different stages 
of the evolution. The solid blue line represents the NIRVANA 
solution, while the red dashed line shows the linear calculation. 

Moving from the top left to top right panel, we can observe 
how the initial pulse begins to broaden. At time t = 0.96 we 
see that the pulse has separated into an in-going and out-going 



bending wave. By time t = 1.46 the pulses have reached the 
boundaries of the domain and have fully separated. It is clear 
that NIRVANA is able to capture the evolution of these low 
amplitude bending waves with a high level of accuracy. 

In lNelson & Papaloizou ( 1999b the numerical SPH and lin- 
ear solutions showed more substantial disagreement at the po- 
sition of the initial pulse (r = 5) at the final stage of the simula- 
tion, with the separation of the two pulses not being captured as 
accurately. This was probably due to the larger numerical diffu- 
sion present in the SPH code, though it should be remembered 
that the SPH simulations were performed using only 20,000 
particles. The more accurate solution obtained with NIRVANA, 
on the other hand, is a clear indication that the numerical diffu- 
sion is smaller. 

We also performed simulations for other values of the disc 
thickness h and maximum initial inclination angle 8 max- The 
number of pressure scale heights used in the meridional direc- 
tion was scaled with the maximum inclination angle to avoid 
the disc interacting with the meridional boundaries. The results 
are shown in Fig|2l and are presented at a stage of the simu- 
lations when the bending waves have propagated all the way 
to the radial boundaries. This corresponds to the final panel of 
FigHI Since the bending waves propagate with half the sound 
speed, the evolution occurs over a longer time in the thinner 
discs (h = 0.01), and over a shorter time in the thicker disc 
models (h = 0.05). As can be seen from Figf2] the agreement 
between the numerical solution and the linear calculation is 
good in all the models presented here. This is particularly the 
case for those runs with lower initial maximum inclinations in 



M.M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 



9 



Fig El (top left and right panels), but the results are still in very 
reasonable agreement for the higher inclination cases in Figj2] 
(lower left and right panels). 

We mention here that for even higher inclination angles 
NIRVANA behaves more diffusively, and the agreement be- 
twee n the linear and non linear solutions gets worse. As found 
by (INelson & Papaloizoull 19991) . when the warps become non- 
linear they generate horizontal motions that are of the order 
of the sound speed, creating shocks because of the symmetric 
form of the initial bending disturbance. The result is an effec- 
tively higher viscosity operating in the disc, which is clearly 
not accounted for in the linear calculations. 

5. Tilted discs in binary star systems 

We now discuss our results for the misaligned, tidally interact- 
ing discs. We set up the disc models according to the procedure 
outlined in Sect. 3, and details of the model parameters can be 
found in Table 1 . The models are characterized by the disc as- 
pect ratio, h, viscosity, a, and initial inclination with respect to 
the binary orbit plane, yp. The disc outer edge, R, was chosen 
such that the disc is already very close to its tidal truncation 
radius, which is about one third of the distance between sec- 
ondary and primary star, D. In some models the initial radius 
of the outer disc edge was varied also. 

In this study, we are interested in examining the evolu- 
tion of discs that fulfil the condition for wave-like warp prop- 
agation, h > a, as well as discs that support diffusive warp 
propagation, h < a. For the models in the wave regime we 
set a = ~h, and for the runs in the diffusive regime we set 
a = 0.1 » h. 

The perturber is evolved on a circular orbit at a distance 
of D — 30, and its mass is increased linearly from zero up to 
Ms = Mp over a time interval of 4 orbits, during which time 
its distance is kept constant. Models 1-5 were performed in 
the precessing frame. If we use Eq. (TTBT l to estimate the preces- 
sion frequency, then we obtain a value of Qf = -4.86 ■ cos(jf) 
[degrees/orbit]. We find, however, that a better fit to test simu- 
lations is given by Q. F = -3.74 • cos(jf) [degrees/orbit], which 
is the value used in the simulations. The reason is that as the 
disc gets tidally truncated it tends to precess slower, as can be 
seen from Eq. (TT3T >. 

The resolution was chosen such there are 5 cells per pres- 
sure scale height in the meridional direction for all models. 
Models 1 and 3 incorporated 15 scale heights in the meridional 
direction. Models 2, 4 and 5 used 22.5 scale heights (to allow 
for a greater degree of twisting in the higher viscosity mod- 
els where warp communication is expected to be less efficient). 
In the radial and azimuthal directions we used N r = 404 and 
N,p = 300 cells, respectively, for each of these models. 

Models 6 and 7 were performed in the non-precessing bi- 
nary frame. Since the disc midplane is inclined with respect to 
the equatorial plane of the computational grid in these simu- 
lations, higher resolutions are necessary to avoid numerical er- 
rors. This is particularly true in the azimuthal direction because 
tilting the disc causes a component of the disc vertical struc- 
ture to point along this direction. For a disc whose midplane 
coincides with the equatorial plane of the computational grid, 



pairwise cancelation of fluxes ensures conservation of angular 
momentum to machine accuracy. However, if the disc midplane 
is inclined with respect to the equatorial plane of the computa- 
tional grid, the accuracy is limited by the advection scheme, 
and higher resolution in the azimuthal direction becomes nec- 
essary to avoid unphysical evolution of precession and inclina- 
tion angles. Hence we used 600 cells in azimuth. In order to 
be able to resolve spiral waves for these extremely thin discs, 
we used 1056 cells in radius. The disc outer edge was chosen 
to be a bit smaller than in the thick disc runs, since highly pro- 
nounced non linear effects at the outer edge were seen in lower 
resolution test simulations (described as Models 6a and 7a in 
Table 1). In effect a strongly perturbed and narrow outer rim of 
gas became partially detached from the main body of the disc 
and developed a very different inclination and precession angle 
evolution in Model 6a, an effect not observed in Model 6 with 
a smaller initial outer radius. This phenomenon is discussed in 
more detail later in the paper. In order to accommodate an in- 
clined disc with an inclination of 10° we used 60 pressure scale 
heights (~ ±17.03 degrees) in the meridional direction, with 5 
cells per pressure scale height. 

We now describe the results of these models, where we 
have ordered our discussion in terms of the disc aspect ratio, 
h. 



5.1. Models 1 and 2 

A column density plot for the disc in Model 1 is displayed in 
the upper panels of Figj3] corresponding to the end stage of 
this simulation. The left panel displays a projection of the disc 
column density onto the (e2, e3) plane of the binary reference 
frame, and the right panel displays a projection onto the binary 
(ei, e^) plane. At this stage the disc has precessed rigidly to 
(3 = -180 degrees, so the disc appears edge on in the (e2, e3) 
plane and face on in the (ei, e3) plane. The precession of 180 
degrees is evident, as the disc angular momentum vector now 
has a negative e2 component, whereas it had a positive e.2 com- 
ponent in its initial state, as can be seen from the transformation 
given by Eq. (|8]l at t — 0. Another result of the interaction with 
the secondary star is the formation of spiral waves, which are 
apparent in Fig|3] 

Fig0]shows the evolution of the precession (J3) and inclina- 
tion (5) angles, respectively, for disc material orbiting at differ- 
ent radii. These have been calculated using Eqs. ( l22l and ( |2"TT i 
respectively, where L has been calculated using Eq. (fT9l , which 
we have integrated over radial shells of width Ar = 1 . Thus the 
red dashed line in Fig|4]corresponds to disc gas between r — 1 
and r — 2, the blue solid line to disc gas between r — 8 and 
r = 9, and the black dotted line to disc gas between r — 4 and 
r = 5. The solid black line shows the precession and inclina- 
tion angles of the entire disc. It is apparent that the lines are al- 
most indistinguishable in Fig0](upper left panel), showing that 
the different disc parts in Model 1 precess at the same uniform 
rate, such that the disc is in a state of rigid body precession 
with al most no twist. This behaviou r is expected from linear 
theory (IPapaloizou & Terquemlll995l) . In Table 1 we show the 
warp propagation timescale t w and the differential precession 
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Model 


h 


a 


7f 


R 


Frame 


Tw 


Tp 


Predicted 


Observed 


label 
















behaviour 


behaviour 


1 


0.05 


0.025 


45 


9 


precessing 


5.43 


48.01 


rigid precession 


rigid precession 


2 


0.05 


0.1 


45 


9 


precessing 


10.87 


47.60 


rigid precession 


rigid precession 


3 


0.03 


0.015 


45 


9 


precessing 


9.05 


43.90 


rigid precession 


rigid precession 


4 


0.03 


0.1 


45 


9 


precessing 


30.19 


47.02 


rigid precession 


rigid precession 


5 


0.03 


0.1 


25 


9 


precessing 


30.19 


36.37 


rigid precession 


rigid precession 


6 


0.01 


0.005 


10 


8 


binary 


22.77 


33.87 


rigid precession 


rigid precession 


7 


0.01 


0.1 


10 


8 


binary 


227.7 


36.25 


disrupted/broken 


rigid precession 


6a 


0.01 


0.005 


10 


10 


binary 


31.83 


21.48 


disrupted/broken 


broken 


7a 


0.01 


0.1 


10 


10 


binary 


318.32 


19.12 


disrupted/broken 


disrupted 



Table 1. Table of runs. The disc is characterized by the aspect ratio, h, viscosity parameter, a, and initial inclination angle yp. 
We also list the initial disc outer radius R and the reference frame in which the simulation was performed. The warp propagation 
times, tw and differential precession times, Tp, for each of the models are also listed and indicate the theoretically expected and 
observed behaviours. 
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Fig. 3. Column density plot at time t — 75.0 for Model 1 (upper panels), and Model 2 (lower panels). The left panels correspond 
to projection onto the (e2, e3) binary plane, and the right panels to projection onto the (ej, 63) plane. 
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Model 1 




Model 2 




Model 1 



47 r 




Model 2 



47 r 




Fig. 4. Evolution of precession (left panels) and inclination angles (right panels) for Models 1 and 2. The blue solid line corre- 
sponds to disc material between r — 8 and r — 9, the red dashed line to disc material between r = 1 and r = 2, and the black 
dotted line to disc material between r = 4 and r = 5. The black solid line represents the entire disc. 



timescale, t>. In the bending wave propagation regime (h > a), 
Tw = 2R/cs, since warps propagate with half the sound speed 
in this regime. In the diffusive regime (h < a) tw — 2a 2 R 2 /v. 
The warp propagation timescale should be compared to the dif- 
ferential precession timescale t> ~ 2/Clp. If tw < ~cp it is 
expected that the disc will precess as a rigid body, while if 
Tw > Tp the disc may become disrupted as a result of differ- 
ential precession. Hence, we see that the rigid precession of 
the disc in Model 1 agrees with expectations. Ou r Model 1 is 
very similar to Model 1 in lLarwood et al. (1996), which also 
shows rigid body precession. The inferred precession period in 
our Model 1 is about 146 orbits, which is a little larger than the 
period given by Eq. (flSb . The difference arises because the disc 
is tidally truncated and shrinks slightly, so its precession rate 
decreases. 

Looking at the upper right panel of Figj4] we observe that 
the inclination for Model 1 changes only slightly during the 
simulation. The oscillation seen in the inclination rates is driven 
by the secondary star. Since the orbital plane is inclined with re- 
spect to the disc midplane, the vertical component of the grav- 
itational force due to the secondary star causes the disc to per- 
form a forced o scillation with a frequency equa l to twice the bi- 
nary frequency. IPapaloizou & Terquerr] ( 1995 ) argued that the 
disc is expected to align with the binary orbit on a time scale 
similar to that over which binary torques cause the total angu- 
lar momentum of the disc to evolve. If the viscously-induced 
outward expansion of the disc is balanced by tides due to the 



companion, then disc-alignmnent should occur on the viscous 
evolution time, T v . For Model 1 we estimate T v = R 2 /v ~ 2175 
orbits. Extrapolating the disc evolution shown in Fig. |4] gives 
an alignment time of ~ 2025 orbits, in good agreement with 

the viscous time scale . 

Lubow & Ogilvie (2000) present images which represent 



disc shapes for a variety of disc parameters, obtained using a 
linear analysis of tilted discs in binary systems. They introduce 
a parameter e, which is a measure of the disc thickness, and is 
roughly equivalent to y[6H/R. Their model a displayed in their 
Figure 7 has H/R a 0.041, a = 0.01 and D/R = 0.3, similar to 
our model 1 . Although a detailed comparison is not possible, 
it appears that there is general agreement between linear the- 
ory and our non-linear simulation, since both show almost no 
distortion of the disc due to twisting or warping. 

Column density plots for Model 2 are displayed in the 
lower panels of Figf3] Warp propagation in this higher viscos- 
ity model (a = 0.1) is expected to be diffusive, and therefore 
less efficient than for Model 1 (as shown by the warp propa- 
gation times listed in Table 1). This run was also evolved until 
the precession angle reached -180 degrees. The high viscosity 
operating in the disc leads to the damping of the spiral waves 
before they can propagate very far, and thus they are not appar- 
ent in the images. 

The precession angles for Model 2 are displayed in Figj4] 
(lower left panel). We observe that this more viscous disc de- 
velops a small differential twist before it attains a state of 
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Fig. 5. Column density plots for Model 3 (upper panels) and Model 4 (middle panels) at time t = 35.0, and Model 5 at time 
t = 55.0. The left panels correspond to projection onto the (ea, 63) plane, the right panel to projection onto the (ei, 63) plane. 



rigid precession. This suggests that the disc needs to develop 
a slightly more distorted shape in order for internal stresses to 
counterbalance the companion-induced differential precession 
when warp communication is less efficient. Nonetheless, rigid 
body precession is expected in this case, and the numerical re- 
sults are in agreement with this expectation. 



Examining the inclination evolution for Model 2 shown in 
the lower right panel of Fig. H] we see it is very similar to 
that observed for Model 1. The viscous time scale in Model 
2 should be approximately 4 times shorter than for Model 1, 
but there is no indication that the inclination evolution is faster 
in this case. We note that the inclination evolution rates do not 
appear to have reached final steady values for either Models 1 
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or 2 by the end of the simulations, such that an accurate assess- 
ment of the longer term evolution rate is not possible. This may 
be because the outer disc radius has not had time to equilibrate 
during the simulations as they have only run for a fraction of the 
global viscous timescale of the disc. We are able to conclude, 
however, that the inclination evolution occurs over timescales 
much longer than the precession time, and appears to be con- 
sistent with evolution on the global viscous time of the disc. 

5.2. Models 3, 4 and 5 

As shown in Table 1, Models 3, 4 and 5 all have h = 0.03. 
Model 3 has a = 0.015, and so warp propagation occurs via 
bending waves, whereas Models 4 and 5 have a = 0.1 such that 
warps propagate diffusively. Table 1 also shows that the warp 
propagation for Models 4 and 5 is expected to be considerably 
less efficient than for Model 3. Models 3 and 4 had inclinations 
of 45°, whereas Model 5 had an inclination of 25°, which in- 
duces a more rapid precession of the disc, and therefore has a 
greater tendency to twist the disc up. 

At the outset of this project it was our intention to run all 
models until they had precessed by 180°. For models 3 and 
4, however, we found that the discs become eccentric on a 
timescale of approximately 60 orbits (16 binary orbits), lead- 
ing to undesirable changes in the disc structure due to our use 
of an open inner boundary condition. Thus we did not evolve 
these models until they had precessed by 180°. We note that the 
timescale for the observed di sc eccentricity growth correspond s 
closely to that obtained by (IKlev. Papaloizou & Ogilvid l2008). 
Model 5 did not suffer from this problem because the lower 
inclination angle induced more rapid precession, such that this 
model could be evolved until it had precessed by 1 80°, but prior 
to growth of significant eccentricity. The study of disc eccen- 
tricity goes beyond the scope of this paper, but the results for 
models 3 and 4 suggest that very long term integrations may 
lead to the formation of inclined and eccentric discs. The fact 
that we did not evolve models 3 and 4 through 180° of preces- 
sion means that the upper and middle panels of Fig|5] display 
images which correspond to 90° of precession. 

It is clear from these panels that the disc in Model 3 
has precessed as a rigid body, with almost no twist apparent. 
Consequently, after precessing through 90° the disc appears 
edge on when projected into the (ei, e3) plane and face on 
when projected into the (ea, e3) plane. The middle right panel 
of Fig. [5] however, shows that the disc in Model 4 has devel- 
oped a small twist due to the inner disc not having precessed 
by 90° at the same moment when the outer disc has. The lower 
panels of Fig. show column density plots for Model 5, which 
was evolved until the outer disc had precessed by 1 80° (assisted 
by the fact that the precession is faster in this case). It is appar- 
ent from these panels that the disc has developed a small twist 
since the inner disc has not precessed a full 180° by the end of 
the simulation. 

The precession angles for Model 3 are plotted in the top left 
panel of Fig. [6] The blue solid line corresponds to disc material 
orbiting between radii r — 8 - 9, the red dashed line radii be- 
tween r = 1 - 2, and the black dotted line to radii between r - 



4 - 5. A black solid line is also plotted showing the precession 
angle integrated over the whole disc. The fact that these lines 
effectively lie on top of each other shows that the disc precesses 
as a rigid body with essentially no twist. Rigid-body precession 
is expected in this case since the warp propagation time is less 
than the differential precession time (see Table 1). 

The inclination angles for Model 3 are plotted in the top 
right panel of Fig. [6] The line styles and colours correspond 
to the same disc annuli as described above. It is clear that the 
inclination evolution in this case is very slow. The global vis- 
cous evolution time for this disc is approximately 10,000 orbits, 
and a linear extrapolation of the inclination evolution shown in 
Fig.|6]gives a time of ~ 13500 orbits for the disc to fully align 
with the binary orbit. Clearly the disc inclination is evolving on 
the viscous time in this case. 

The evolution of the precession angles for Model 4 are 
shown in the left middle panel of Fig [6] The line styles and 
colours correspond to disc annuli as described above. We see 
that this disc develops a well-defined twist which is set up 
within the first 5-10 orbits due to the inner and outer parts of 
the disc undergoing differential precession, with the outer disc 
precessing faster than the inner disc. The magnitude of the twist 
is approximately 12°. Thereafter the disc is seen to precess as 
a rigid body, maintaining the same twisted shape. We note that 
this disc is expected to achieve rigid body precession from con- 
sideration of the warp propagation and differential precession 
times (see Table 1). 

The inclination angles for Model 4 are shown in the right 
middle panel of Fig. [6] We see that the disc in this case has de- 
veloped a small warp, with the inner and outer discs having a 
difference of inclination of about 1°. We also observe that the 
inclination is evolving more rapidly than is the case for Model 
3. The global viscous evolution time for this disc is approxi- 
mately 1500 orbits, and extrapolating the inclination evolution 
observed in Fig. [6] gives an estimated time for alignment of 
~ 1000 orbits, close to the expected value. 

It is interesting at this point to compare our results for 



Models 3 and 4 with Model 9 presented by lLarwood et al 



(119961) . which had comparable parameters (h = 0.03, D/R = 3, 
45 degree inclination). Although it is difficult to know the pre- 
cise value of the a viscosity operating in the SPH simulations, 
it is likely that the effective a value lies between the values 
used in our Model 3 (a = 0.015) and Model 4 (a = 0.1). 
Nonetheless, we find qualitatively different behaviour, with our 
models quickly achieving a state of rigid body precession, with 
a small twi st and warp which va ry smoothly across the disc. 
Model 9 of lLarwood et al.1 dl996) shows significant disruption 
of the disc, which effective breaks discontinuously into two dis- 
tincts parts which become mutually inclined by approximately 
25°. The origin of this discrepancy is unclear, but one possi- 
bility is that SPH simulations using 20,000 particles are only 
marginally able to resolve the vertical structure of a disc with 
h = 0.03. This may have the effect of reducing the efficacy of 
warp propagation below that which is expected from consider- 
ation of the disc physical parameters. 

We now turn to Model 5. The precession angles for this 
run are displayed in the lower left panel of Fig|6] Given that 
the precession frequency, given by Eq. (fTBT l, is proportional 
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Fig. 6. Evolution of precession (left panels) and inclination angles (right panels). The blue solid line corresponds to disc material 
between r — 8 and r — 9, the red dashed line to disc material between r — 1 and r = 2 and the black dotted line to disc material 
between r — 4 and r = 5. The solid black line represents the entire disc. In the lower left panel the free particle precession 
rates for the inner (dashed-dotted line) and outer disc (dashed-triple-dotted line) are also depicted for Model 5. These have been 
calculated using Eq. 



to cos (jf), the disc precesses faster in this simulation than in 
Models 1-4, and this faster precession can be observed in 
Fig. [6] The lower left panel of this figure shows the preces- 
sion angles as a function of time for disc annuli located be- 
tween r = 1 — 2, r = 4 - 5, and r = 8 — 9, with the linestyles 
and colours being the same as described above. The larger 
precession frequency induces greater differential precession in 
the disc, and the consequence of this is that Model 5 shows 
a slightly larger twist (15°) than Model 4 (12°). Also plotted 
in this panel are precession angles that would be observed if 
the inner (black dashed-dotted line) and outer (black dashed- 
triple-dotted line) disc annuli precessed at their free particle 



rates given by (Papaloizou & Terq uemll995 ): 
3 GM S 



-(3cos z 5- 1). 



(38) 



m K D 3 

As can be seen from the figure, the disc parts tend initially 
towards their free particle rates, such that the outer disc pre- 
cesses faster than the inner disc. Information about this differ- 
ential precession is communicated across the disc, and internal 
stresses are established which cause the inner disc precession to 
speed up and the outer disc precession to slow down, such that 
precession at a uniform rate occurs after approximately 10 or- 
bits. The warp communication and differential precession times 
given in Table 1 for this model predict rigid-body precession, as 
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Fig. 7. Precession rates due to the contributions discussed in Sect. 3.6 for Model 5. The blue solid line corresponds to disc material 
between r — 8 and r — 9, the red dashed line to disc material between r = 1 and r = 2 and the black dotted line to disc material 
between r — 4 and r = 5. These rates have been calculated using Eqs. ([34-b . In the lower right panel the total rate is displayed, 
which has beeen calculated using Eq. (l33l . where we subtracted the constant mean precession rate Slf of the precessing frame. 
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Fig. 8. Inclination rates due to the contributions discussed in Sect. 3.6 for Model 5. The blue solid line corresponds to disc 
material between r = 8 and r = 9, the red dashed line to disc material between r = 1 and r — 2 and the black dotted line to disc 
material between r — 4 and r = 5. These rates have been calculated using Eqs. (1341 , 



observed. The emergence of internal stresses as the disc twists 
up is illustrated in FigfT] where the different physical contribu- 
tions to the precession rate are displayed for Model 5. These 
have been calculated according to the procedure outlined in 
Sect. 3.6. Here the different linestyles correspond to the same 
disc annuli described above (solid blue line: r — 8 - 9; dotted 
black line: r = 4 - 5; dashed red line: r = 1 — 2). From the lower 



left panel of FigJT] we can observe that the gravitational inter- 
action with the secondary causes a negative precession rate for 
the outer annulus, and a positive precession rate for the inner 
annulus, and these rates correspond to the free particle preces- 
sion rates. Since the frame in which these rates were measured 
precesses in a retrograde sense with a mean frequency Q. F < 0, 
this means that the inner annulus would lag behind while the 
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Fig. 9. Time averaged inclination change rates as function of 
radius for Model 5. Red dashed-dotted line: radial flux of angu- 
lar momentum; blue dotted line: gravitational interaction with 
the companion star; green dashed line: viscous friction between 
adjacent radial disc shells; black dashed-triple-dotted line: total 
inclination change rate. 



outer annulus would lead the precession of the reference frame, 
if the companion's gravity was the only contributor to the to- 
tal precession rate. In fact, this is what is observed until about 
t — 10 orbits, during which time the disc precesses differen- 
tially (see Fig|6] lower left panel). As the disc starts to develop 
a twist due to differential precession, misalignment of disc mid- 
planes as a function of radius causes pressure forces to generate 
horizontal shear motions that drive the propagation of bending 
disturbances, and these communicate information about pre- 
cession angle changes across the disc. The resulting internal 
stresses are clearly proportional to the degree of twisting, al- 
lowing the disc to achieve a state of rigid body precession once 
a sufficiently twisted disc has been established. At this point the 
disc structure becomes quasi-static in the precessing frame as 
each disc annulus precesses at the same rate. This sequence of 
events is demonstrated by the upper left panel of Fig|7] which 
shows the contribution to the precession rate from the pressure- 
induced radial flux of angular momentum. After t — 10 orbits, 
the contribution to the precession rate from the radial flux al- 
most completely compensates the effect of gravity, such that 
the total precession rate (relative to the precession of the ref- 
erence frame) is zero for all disc annuli, as can be seen in the 
lower right panel of Fig|7] which shows the sum of the contri- 
butions to the precession rate. This demonstrates that the disc 
manages to redistribute the angular momentum such that rigid 
body precession is achieved. Once again we note that in this 
figure we have subtracted the mean precession rate Of of the 
precessing frame. 

The effect of viscosity as calculated in Sect. 3.6 is depicted 
in the upper right panel of FigJT] We observe that the effect 
of viscous friction between differentially precessing radially 
adjacent disc shells causes the disc to precess more rigidly. 
However, the dominant effect of the viscosity is the damping of 
the horizontal shear motions induced by the twist. Since these 



shear motions are responsible for driving the bending distur- 
bances, the viscosity will lead to weakened communication, 
and the redistribution of angular momentum due to radial ad- 
vective fluxes shown in the upper left panel of Fig|7]will be sup- 
pressed. This effect is dominant and opposite to the one seen in 
the upper right panel of Fig f7] 

The inclination angles for Model 5 are shown in the lower 
right panel of Fig. [6] It is clear that the disc develops a small 
warp, with the inner disc having ~ 1° higher inclination than 
the outer disc. The rate of global inclination change is similar 
to that observed in Model 4, suggesting alignment of the disc 
in the global viscous evolution timescale. 

We plot the different physical contributions to the rate of 
inclination change for Model 5 in FigJS] In Fig(9]we also dis- 
play the time averaged contributions to the rate of inclination 
change as a function of radius, where the time averaging was 
performed over the full length of the simulation. As may be ob- 
served in the lower left panel of Fig[8] and even more clearly 
in Fig(9](blue dotted line), the companion's gravity induces a 
globally negative net inclination change, leading to coplanarity 
of the entire disc on the viscous time scale, in agreement with 
the low er right panel of Fig|6] This is consistent with the find- 
ings of Papaloizou & Terquem ( 1995) that the zero-frequency 
(secular) gravitational interaction will lead to coplanarity. We 
note that because of the time averaging, this is the dominant 
gravitational contribution that we pick up in Fig |9] (blue dotted 
line). 

More interestingly, we can observe from the upper left 
panel of Fig[8] and in Fig(9] (red dashed-dotted line), that the 
contribution from the radial fluxes not only counteracts the ef- 
fect of gravity, it actually overshoots slightly, causing the inner 
disc to become more inclined than the outer disc. This change 
of inclination is seen in the upper right panel of Fig|6] Thus the 
communication of bending disturbances on the one hand leads 
to rigid precession of the disc, and on the other hand causes the 
disc to become slightly warped. Viscosity tends to reduce the 
amount of warping, such that it acts to flatten the differential 
inclination profile across the disc radius, as can be seen from 
the green dashed line in Fig (9] Thus the total average inclina- 
tion rate in Fig|9] (dashed-triple-dotted line) is negative for the 
entire disc, while the magnitude is slightly bigger for the outer 
disc, with the consequence that the disc has developed a slight 
warp and tends to become coplanar with the orbital plane on 
the viscous timescale, as seen in the lower right panel of Fig|6] 



5.3. Models 6 and 7 

We now discuss Models 6 and 7, for which the discs have as- 
pect ratio h - 0.01, an inclination angle of 10°, and were per- 
formed in a non precessing frame. Model 6 has a = 0.005, 
so warps propagate via bending waves in this case, and Model 
7 had a = 0.1, so warps propagate diffusively in this model. 
We see from Table 1 that differential precession is expected 
to disrupt the discs in Model 7, as the warp propagation time 
is longer than the differential precession time. Rigid-body pre- 
cession is expected for Model 6. 
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Fig. 10. Column density plots for Model 6. The left-hand panels are projections onto the £3) plane, and the right-hand panels 
are projections onto the (ei, 63) plane. The disc structure is displayed at the different times indicated. 
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Fig. 11. Evolution of the precession (left panels) and inclination angles (right panels). The blue solid line corresponds to disc 
material between r — 7 and r = 8, the red dashed line to disc material between r — 1 and r = 2, and the black dotted line to disc 
material between r — 4 and r = 5. The solid black line represents the entire disc. 



Column density plots showing the disc in Model 6 are pre- 
sented in Fig [10] for different times during the evolution. The 
top panels show the initial state of the disc. The disc appears 
edge on in the (e2, e3) plane and face on in the (ei, e3) plane. 
As we move down to the next panels the outer edge of the disc 
has precessed by 45 degrees at time t = 10.4, as can be seen 
in Fig[TT](upper left panel). The inner region of the disc, how- 
ever, has only precessed by about 12 degrees at this stage and 
so appears almost edge-on in the left panel. In the third panels 
down the outer disc has precessed by 90 degrees at t = 25.7 
and appears edge on in the (ei, e3) plane and face on in the (e2, 
e3) plane. The inner disc has precessed by approximately 70 
degrees by this time, and so does not appear edge-on in either 
projections. The bottom panels show the disc structure at the 
final output of the simulation which was halted at t — 44 orbits, 
by which time the outer disc had precessed by approximately 
160°. 

The precession angles for Model 6 are shown in the upper 
left panel of Fig. [TTJ The blue solid line corresponds to disc 
material orbiting in an annulus between r =7-8, the black 
dotted line corresponds to the annulus between r = 4 — 5, and 
the red dashed line corresponds to the annulus between r = 1 
- 2. As with Models 4 and 5, we see that the disc inner and 
outer annuli begin to precess at the local free-particle rates, and 
the disc develops a significant twist which reaches a maximum 
amplitude of approximately 30°. As the disc becomes increas- 
ingly twisted, however, internal stresses are established which 



cause the inner disc precession rate to increase and the outer 
disc precession rate to decrease. Eventually the disc reaches a 
state of rigid body precession after a time of about 10 orbits. 
After this time we see that there is a long term readjustment of 
the degree of twist such that by t — 44 orbits the twist angle has 
reduced from ~ 30° to ~ 20°. 

In Fig. [12] we display the time integrated contributions to 
the precession rate as a function of radius for Model 6, follow- 
ing the procedure described in Sect. 3.6. Note that the average 
(negative) precession rate of the disc has been subtracted in 
this figure, such that a positive value corresponds to preces- 
sion that will lag that of the overall disc, and a negative value 
corresponds to precession that leads. The contribution from the 
companion's gravity is shown by the blue dotted line, and this 
is very similar to the free-particle precession rate which is plot- 
ted using the black solid line. The red dashed-dotted line shows 
the contribution from the pressure-induced radial flux of an- 
gular momentum, and this is seen to almost exactly counter- 
balance the effect of gravity, showing that it is the pressure- 
induced radial motion in the disc which is responsible for halt- 
ing the differential precession and inducing rigid body preces- 
sion. The contribution from the disc viscosity is shown by the 
green dashed line, and is seen to be negligible. 

The inclination angles for Model 6 are shown in the top 
right panel of Fig. Q~T] where the line styles and colours cor- 
respond to the disc annuli described above. As with Model 3, 
we see that the disc develops a warp such that the inner disc 
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Fig. 12. Average precession rates as function of radius for 
Model 6. Red dashed-dotted line: radial flux of angular mo- 
mentum; blue dotted line: gravitational interaction with the 
companion star; green dashed line: viscous friction between 
adjacent radial disc shells; black dashed-triple-dotted line: to- 
tal inclination change rate. The black solid line represents the 
precession rate of free particles. 




Fig. 13. Average inclination change rates as function of radius 
for Model 6. Red-dashed-dotted line: radial flux of angular 
momentum; blue dotted line: gravitational interaction with the 
companion star; green dashed line: viscous friction between ad- 
jacent radial disc shells black dashed-triple-dotted line: total 
inclination change rate. 



has a larger inclination to the binary plane than the outer disc. 
Indeed, the inclination of the inner disc is seen to increase ini- 
tially, before decreasing after approximately 20 orbits as the 
disc as a whole aligns with the binary orbit plane. As discussed 
previously, disc alignment is expected to occur on the viscous 
time scale, which for Model 6 is approximately 3 x 10 5 orbits. 
It is obvious from Fig.Qj] however, that the disc in our simula- 
tion is not aligning on such a long time scale, but will actually 
align within ~ 650 orbits. Possible reasons for this enhanced 
alignment rate are discussed below. 



Fig- H3 shows the time integrated contributions to the rate 
of inclination evolution as a function of radius for Model 6. 
We see that the companion's gravity causes the inclination to 
decrease in the disc inner and outer regions, and in the inner re- 
gions the radial flux of angular momentum opposes this and in- 
creases the inclination during the simulation. Indeed the radial 
flux contribution causes the inclination angle to grow there, as 
observed in the lower right panel of Fig. QT| The outer parts of 
the disc experience negative contributions from both the com- 
panion's gravity and the pressure-induced radial flux, driving 
alignment of the disc on the relatively short timescale described 
above. It is clear that the disc viscosity plays a negligible role 
in driving the inclination evolution (green dashed line). 

We may compare our model 6 with model c in Figure 7 of 
Lubow & Ogilviel d2000l) . which had H/R =* 0.012, a = 0.01 



and D/R _ 0.3. Detailed compari son between the results pre- 
sented in iLubow & Ogilv ie (2000) and our model 6 is not pos- 
sible, but we note that there is general agreement in the pre- 
diction that a disc with H/R 0.01 and a < H/R will show 
significant distortion due to twisting and warping. A detailed 
comparison between our non linear simulation results and the 
predictions of linear theory will be presented in a future publi- 
cation. 

We now discuss the results from Model 7. In Fig[l4]we dis- 
play column density plots for this model. The top panels show 
the initial state. The second panels down show the disc at time 
t = 20.5. At this stage the outer disc edge has precessed by 90 
degrees and appears edge on in the (ei, e3) plane, and face on 
in the (e2, e3) plane. By contrast, the inner disc has not pre- 
cessed very much and still appears almost edge on in the {ez, 
63) plane and face on in the (ei, e3) plane. At time t — 44 orbits 
the outer disc has precessed by 180 degrees and appears edge 
on in the (e2, e3) plane and face on in the (ei, e3) plane. The 
inner disc, however, has only precessed by about 70 degrees 
by this time, and so does not appear edge on in either projec- 
tion. The bottom panels show the final stage of the simulation 
when the outer disc has precessed by about 220 degrees. It is 
clear from these images that the disc in Model 7 has adopted a 
highly twisted shape, but one in which the precession angle for 
different disc annuli varies very smoothly across the disc. 

The precession angles for Model 7 are displayed in the 
lower left panel of FigQj] We observe that the disc is in a state 
of differential precession at the end of the simulation, and the 
twist between inner and outer disc annuli is approximately 110 
degrees at this stage. Although this twist is large, the actual rate 
of differential precession at the end of the simulation is very 
small, as the precession rates of the inner and outer disc have 
converged during the run. Indeed, extrapolation of the preces- 
sion angles in Fig. fTTI indicates that rigid body precession will 
be achieved at a time equal to t — 56.6 orbits, at which point the 
twist amplitude will be 1 1 1.7°. This suggests that even though 
the warp propagation time is much smaller than the differential 
precession time for Model 7 (see Table 1), internal stresses can 
be set up within the disc that cause uniform precession once 
the disc becomes highly twisted. Adopting a definition of disc 
disruption which requires the twist to become greater than 180 
degrees, we suggest that the disc in Model 7 will not be dis- 




Fig. 14. Column density plots for Model 7. The left-hand panels are projections onto the £3) plane, and the right-hand panels 
are projections onto the (ei, 63) plane. The disc structure is displayed at the different times indicated. 
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Fig. 15. Diagram illustrating the alignment torques which arise 
in highly twisted discs, as discussed in the text. 

rupted, but will instead adopt a highly twisted but uniformly 
precessing configuration. 

The inclination angles for Model 7 are shown in the lower 
right panel of Fig. QT| where the lines styles correspond to the 
disc annuli described above. We can see that the disc in this 
case develops a small warp with the difference between the in- 
clinations of the inner and outer disc being ~ 1°. We also note 
that it appears that the integrated tilt of the whole disc is ac- 
tually smaller than for any of the individual disc annuli, but 
this is an effect of averging over a significantly twisted disc. As 
observed for Model 6, we see that the inclination of the disc 
evolves rather quickly. The global viscous evolution time for 
Model 7 is ~ 15000 orbits, whereas extrapolation of the in- 
clination angles for either the inner or outer disc suggests that 
they will completely align after ~ 100 orbits. 

What is the explanation for this rapid alignment observed 
for Models 6 and 7 ? One possibility that we have explored is 
that numerical diffusion may cause a tilted disc to align with 
the equatorial plane of the computational grid, since the advec- 
tion routine does not conserve total angular momentum when 
the disc azimuthal velocity is not directed along the azimuthal 
direction of the computational grid. Lower resolution test cal- 
culations performed for a tilted disc without a companion indi- 
cate that this effect can cause eventual disc alignment, but at a 
rate which is at most 7 times smaller than observed in Models 
6 and 7. This suggests that numerical diffusion is not the cause 
of the rapid alignment. 

Another possibilty arises when we consider the evolution 
of a rigidly precessing, tilted disc which is in a precessing ref- 
erence frame whose precession frequency is equal to that of the 
disc. Consider first the situation in which the angular momen- 
tum vectors for all disc annuli lie in the y-z plane of the binary 
frame initially (i.e. the z direction of the precessing frame lies 
initially in the y-z plane of the binary frame, pointing along 
the positive y and z directions). A torque exerted on disc annuli 
in the y direction, (a 'y-torque), will cause the disc inclination 
to change. A torque exerted in the x direction (an 'x-torque'), 
will cause precession at a rate which differs from the preces- 
sion rate of the frame. A positive x-torque causes precession 
to occur at a rate which is faster than the frame, and a nega- 
tive x-torque causes slower precession. This general picture is 



expressed mathematically by Eqs. d34l >. and is shown diagra- 
matically in the left hand image of Fig. Q3] For a disc which 
precesses uniformly without any twist, then the companion's 
gravity will exert positive x-torques on the outer disc annuli, 
and negative x-torques on the inner disc annuli, as it tries to 
cause the disc to precess differentially. Internal stresses, how- 
ever, will exactly counterbalance these x-torques to maintain 
uniform precession. In this simple scenario there is no torque 
exerted on the disc in the y direction to modify its inclination. 

Consider now the case of a highly twisted disc which is be- 
ing forced to precess uniformly, and whose annuli all have the 
same inclination. Gravity will again try to induce differential 
precession, but now the torque exerted on disc annuli will have 
both an x and a y component due to the twist. Resolving the 
vectors associated with the internal stresses required to enforce 
uniform precession shows that their combined x-torques oper- 
ate in opposite directions such that they cancel when integrated 
over the disc. The combined y-torques, however, operate in the 
same direction, implying that the internal torques generate a net 
y-torque on the disc. Conservation of angular momentum ob- 
viously prohibits this from arising, suggesting that a y-torque 
must operate, whose effect is to change the disc inclination. 
This is shown diagramatically in the right panel of Fig. Q3] 
Given that the internal torques are operating on the precession 
timescale in order to maintain uniform precession, this argu- 
ment suggests that a highly twisted disc may align on the pre- 
cession time, as is observed for Models 6 and 7. Interestingly, 
the resultant y-torque that arises in this picture is proportional 
to sin (A/3), where A/3 is the twist amplitude, and the rate of in- 
clination evolution observed in Models 6 and 7 is in agreement 
with this scaling. 

If the above argument is correct, then it suggests that discs 
with modest levels of twist will align on the global viscous 
timescale of the disc, but highly twisted discs may align on 
the much shorter precession timescale. This obviously has sig- 
nificant implications for astrophysical systems which may be 
expected to harbour highly twisted discs, such as the X-ray bi- 
naries Her X-l and SS433. 

5.3.1 . Models 6a and 7a - broken or disrupted discs ? 

We now briefly discuss the results of two lower resolution mod- 
els which were run during an early stage of this project: Models 
6a and 7a. These were identical to Models 6 and 7 except that 
the number of grid cells in the radial direction was N r = 300, 
and the disc outer radii were located at R = 10 instead of R = 8. 
Extending the disc radius outward has the potential effect of in- 
creasing the rate of differential precession across the disc, pro- 
vided that tidal truncation does not cause the disc to shrink back 
down to R — 8. 

Column density plots for Model 6a are shown in the top 
two panels of Fig. [16] corresponding to an evolution time of 
t = 22.8 orbits. We can see that there has been some tidal trun- 
cation of this disc, and the outer disc edge is significantly dis- 
torted in the images. This arises because a narrow outer rim, 
running between radii r ~ 7 and r ~ 10, has effectively de- 
tached itself from the main body of the disc, and has evolved 
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Fig. 16. Column density plots of the two lower resolution runs Models 6a and 7a. The upper panels show two different projections 
for Model 6a at time t = 22.8 The lower panels show two different projections for Model 7a at time t - 27.0. 



to have significantly different precession and inclination angles 
compared to the main body of the disc. 

The inclination and precession angles for Model 6a are dis- 
played in the upper left and right panels of Fig. [17] respectively. 
The precession angle of the annulus which lies between radii 
r = 1 — 2 is shown by the thick red solid line, and the remain- 
ing red dotted lines show the precession angles for disc annuli 
of unit width out to r = 7. The solid blue line corresponds to 
the disc annulus lying between r — 9 - 10, and there are two ad- 
ditional blue dotted lines for annuli lying between r — 7 - 8 and 
r — 8 - 9. It is clear that the disc has separated into two parts 
discontinously at a radius between r — 7 and r = 8, with the 
detached outer rim having precessed significantly faster than 
the inner disc. The inclination angles plotted in the upper right 
panel show that the disc becomes significantly warped, with 



the outer rim tending to align with the binary and the inner 
disc actually increasing its inclination. Although the change in 
inclination appears to occur fairly smoothly across the disc, it 
should be noted that the disc has developed a fairly large warp 
of 5°. 

The behaviour of Model 6a is re miniscent of the broken 
disc presented by lLarwood et al. (1996) (their broken Model 9 
had h - 0.03, inclination of 45°, D/R - 3) except that the dis- 
continuous break in Model 6a occurs near to the edge of the 
disc. It would appear that Model 6a breaks because the larger 
disc radius induces a larger rate of differential precession, and 
also possibly because warp propagation is retarded in the outer 
disc regions which are subject to non linear density structures 
and shocks because of strong interaction with the binary com- 
panion. 
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Fig. 17. Evolution of precession (left panels) and inclination angles (right panels) for Models 6a and 7a. The blue solid lines 
correspond to disc annuli of unit radius lying between r — 7 and r — 10. The solid blue line denotes the annulus lying between 
r — 9 and r = 10. The red lines correspond to disc annuli of unit width lying between r = 1 and = 7. The solid red line denotes 
the annulus lying between r = 1 and r = 2. 



Column density plots for Model 7a are shown in the lower 
two panels of Fig. [16] corresponding to a run time of t = 27 
orbits. Here we can see that the disc has become very highly 
twisted (the twist amplitude in the figure is ~ 150°), but 
does not show any of the discontinuous behaviour observed 
in Model 6a. Instead the twist varies very smoothly across the 
disc. The precession and inclination angles for Model 7a are 
shown in the lower left and right panels of Fig.[T7J respectively. 
The solid red line shows the precession angle of the annulus lo- 
cated between r = 1 — 2, and the solid blue line corresponds to 
the annulus between r — 9 - 10. By the end of the simulation 
the twist amplitude is ~ 170°, and the precession rates of inner 
and outer disc remain very different such that this disc will dif- 
ferentially precess through 180°. So, by our working definition 
of disc disruption, this disc will be disrupted because of differ- 
ential precession. The inclination angles show that the disc also 
becomes significantly warped, with the outer disc aligning with 
the binary orbit more rapidly than the inner disc. 

To summarise: Model 6a and 7a both show examples of 
discs which break or are disrupted because of differential pre- 
cession. In the case where warps propagate via bending waves, 
the disc breaks discontinuously into two independently pre- 
cessing parts, but which themselves remain as coherent struc- 
tures which precess as rigid bodies. In the case where warps 
propagate diffusively, viscosity is able to smooth out any dis- 
continuous behaviour, and instead the disc becomes disrupted 
because the twist exceeds 180 degrees. In all probability the 
disc annuli in Model 7a will continue to precess differentially in 
the longer term, substantially destroying any coherent disc-like 
structure. The larger disc radius in this case, when compared 
with Model 7, causes the differential precession rate across the 
disc to be larger, and thus explains why Model 7a shows dis- 



ruption, whereas Model 7 results in a highly twisted disc which 
is able to attain a state of rigid-body precession. 

6. Conclusions 

We have presented non linear simulations of accretion discs 
in close binary systems in which the binary midplane is mis- 
aligned from the binary orbital plane. Previous work on this 
problem, which employed line ar theory, semi-analytic tech 



niques and SPH simulations (Papaloizou & Terquem 1995, 
Larwood et al. 1 19961: iLubow & Ogilvid 12000b . suggests that 
under suitable conditions the disc should become mildly 
warped and precess as a rigid-body around the angular momen- 
tum vector of the binary system. The main aim of the present 
study is to examine in detail the structure and evolution of 
discs in misaligned binary systems, subject to variation of the 
important physical parameters. We used the grid-based code 
NIRVANA to perform the numerical simulations, and in prin- 
ciple this should have the advantage over previous SPH simu- 
lation studies that the disc viscosity is a well defined and con- 
trollable para meter. 

Following lNelson & Papaloizou d 19991) . the propagation of 
linear bending waves was used as a means to calibrate the code 
and test its ability to propagate warps. Direct comparison with 
linear calculations show that the code can propagate bending 
waves with a high degree of accuracy for a range of disc models 
and warp perturbation amplitudes. 

In our study of discs in misaligned binary systems, a num- 
ber of different models were considered. Discs with aspect ra- 
tios h = 0.05, 0.03 and 0.01 were studied, and for each of 
these values models were computed with viscosity parameters 
a = h/2 and a = 0.1. The smaller value of a allows bending 
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waves to propagate, whereas the larger value causes warps to 
evolve diffusively. 

It is expected that discs for which the warp propagation 
time, Tw, is shorter than the differential precession time, Tp, 
will undergo rigid-body precession. Discs which do not satisfy 
this criterion may be disrupted by strong differential preces- 
sion. All our models with h = 0.05 and 0.03 have Tw < Tp, and 
the simulations show that each of these discs attains a state of 
rigid-body precession. Discs with lower viscosity maintained a 
precessing structure with almost no discernible twist and warp, 
due to the highly efficient warp propagation associated with 
bending waves. Discs with larger viscosity propagate warps 
less efficiently, and these models underwent a short period of 
differential precession, setting up a twist in the disc, prior to 
attaining a state of rigid-body precession. Our analysis shows 
that internal stresses are established in the disc as it becomes 
twisted, which transport angular momentum across its radius 
and cause it to precess rigidly. Models in which warp propaga- 
tion is the least efficient develop the largest twist, since the in- 
ternal stresses which hold the disc together against differential 
precession are proportional to the twist amplitude. Examination 
of the simulations shows that the dominant contributor to these 
internal stresses are pressure-induced radial angular momen- 
tum fluxes which are driven by local misalignment of disc mid- 
planes. When a — 0.1, the final amplitude of the twist was just 
a few degrees for the h - 0.05 disc, and for the h - 0.03 disc it 
was between 12-15 degrees, depending on the binary inclina- 
tion. The low viscosity models showed essentially no twist. 

In addition to becoming twisted, a disc may also become 
warped. We find that for all models with h = 0.05 and h = 0.03 
the degree of warping is very modest, with differences in in- 
clination across all disc annuli being less than 1 degree. It is 
expected that these discs will also align with the binary or- 
bital plane on the global viscous evolution time of the disc 
(Pap aloizou & Terquem 1995 : Larwood et alJ [l 996). and all 
simulations with h = 0.05 and 0.03 are fully consistent with 
this expectation. 

We also considered a number of thin disc models with 
h = 0.01, where we not only varied the viscosity but also the 
outer radius (taking values of either R = 8 or R = 10). The 
low viscosity model with R — 8 developed a significant twist 
of between 20 - 30° before achieving a state of solid body pre- 
cession, as expected since tw < Tp for this model. A model 
run with a = 0. 1 and R = 8 was predicted to be disrupted due 
to differential precession since tw > Tp in this case. During 
the simulation it developed a very large twist ~ 1 10°, but as the 
disc become increasingly twisted the rate of differential preces- 
sion decreased. At the end of the run the disc was experiencing 
very slow differential precession, and extrapolation of the rate 
of twisting indicates that this disc will achieve rigid-body pre- 
cession with a twist of ~ 112°. This is the first indication in 
our results that a disc which is predicted to be disrupted can 
nonetheless generate internal stresses to enforce rigid preces- 
sion when the degree of twist becomes large. 

Models with h = 0.01 and R — 10 showed different be- 
haviour, largely as a result of the larger precession rate which 
is induced for a larger disc, and suggest that the models run 
with R — 8 had outer radii slightly smaller than the natural 



tidal truncation radius. The disc with a = h/2, which supports 
bending waves, was observed to undergo modest tidal trunca- 
tion and to break into two distinct disc parts which precessed 
independently of one another, an outco me which is similar to 
one obtained using an SPH simulation bv lLarwood et al.l (fl996) 
but for a disc with somewhat different parameters. A narrow 
rim at the edge of the disc, running between radii r = 7 - 10, 
detached from the interior disc and precessed independently at 
a faster rate. The inner disc part precessed at a rate similar to 
that observed in the R = 8 case. Interestingly we estimate that 
Tw > Tp for this larger disc model, such that the breaking of 
the disc is indeed expected. 

The model with a = 0.1 propagates warps diffusively, and 
was observed to differentially precess through 180 degrees, but 
maintained a very smooth twist profile throughout the simula- 
tion. This disc is clearly disrupted through strong differential 
precession, and the internal stresses are insufficient for the disc 
to achieve rigid-body precession. Over longer times this disc 
should become completely twisted up by differential preces- 
sion such that it loses all semblance of a disc-like structure. 

The results for the R = 10 discs are interesting, as they 
indicate two different modes by which a disc may break or be 
disrupted. In the bending wave regime, our results suggest that 
disc disruption may occur via a discontinuous change in the 
precession rate at a particular radius, such that the disc breaks 
into two distinct pieces which precess independently. Within 
each separate disc-piece warp communication allows for rigid- 
body precession. The detachment of the disc outer rim may in 
part be influenced by the strong tidal interaction there, since 
non linear features and shocks may have the effect of reducing 
the efficacy of warp communication. In the diffusion regime, it 
appears that viscosity is able to smooth out any discontinuities, 
and the disc is disrupted through global differential precession 
of the disc annuli, with a smooth twist profile being maintained 
across the disc radius. 

A significant result obtained for the h = 0.01 models 
is that these highly twisted discs align with the binary orbit 
plane much more quickly than the thicker discs, which align 
on the viscous timescale. Although we find that numerical dif- 
fusion can cause an inclined disc to align with the equatorial 
plane of the computational domain, this occurs on much longer 
timescales than are observed. We tentatively suggest that align- 
ment of highly twisted discs occurs on the precession time 
rather than the viscous time, causing rapid alignment of very 
thin discs. 

We have not considered different binary mass ratios, instead 
we restricted our study to the case of an equal mass binary. 
As Eqs. ( 1381 ) and (TT~5T > indicate, the free and mean precession 
rates scale linearly with the companion mass for a disc of fixed 
outer radius. Hence we would expect the differential precession 
timescale Tp to be increased for less massive binary compan- 
ions, and consequently the differential twist observed in our 
models should be reduced in such systems. 

Our work has a number of astrophysical implications. Most 
T Tauri stars are members of binary or multiple systems, and in 
a number of these it is believed that the disc and orbit planes are 
misaligned. One particular system for which there is a resolved 
image of a disc which is misaligned from the binary plane is 
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HK Tau dStapelfeldt et al.lll998l) . In this system, the disc radius 
is estimated to be 105 AU, and the projected binary separation 
is approximately 3 times larger than this, suggesting that the 
disc may be tidally truncated and undergoing strong interac- 
tion with the inclined companion star. The images do not show 
any signs that the disc is warped, which is consistent with our 
h = 0.05 models, whose aspect ratio is probably appropriate 
for T Tauri discs. Evidence for misaligned young binaries is 
also provid ed by observations of precessing jets in star form- 
ing regio ns Terquemet al. (1999). In one particular example, 
Ballv & Devinel d 19941) suggest that the jet which appears to 
be launched by the young stellar object HH43* in the L1641 
molecular cloud in Orion precesses with a period of about 10 4 
years. Applying models to a disc with outer edge of 50 AU sur- 
rounding a solar type star, the rigidly precessing disc models 
presented in this work precess with a period of between 3.5 - 
5 x 10 4 years, depending on the inclination angle considered. 
Thus the observed precession is consistent with that driven by 
a binary companion, with parameters close to those adopted 
in this work. The long alignment timescales that we find for 
h = 0.05 discs suggests that a T Tauri disc will remain mis- 
aligned throughout its lifetime, such that jets launched from 
the disc inner regions will continue to precess for the duration 
of the time when jet launching is active. 

The thin h = 0.01 disc models that we have computed are 
most likely relevant to discs in X- ray binaries. The X-ray bi 



nary sy stems SS433 Margon (1984) and Her X- 1 iBoynton et al 



; 1980) are two examples of systems thought to contain pre- 
cessing discs, with the evidence provided by the precessing 
jet of SS433 being particularly compelling. It is likely that 
the transfer of matter between the donor star and compact ob- 
ject in these systems will cause disc material to lie in the bi- 
nary orbit plane initially, and a tilting mechanism is required 
to move the dis c out of this plane, such a s the Bardeen- 



Petterson effect (Bard een & Petterson 1 119751) . a misaligned 



dipole magnetic field associa ted with the c entral neutron star 
dTerquem & Papaloizoull2000 ), a disc wind dSchandl & Mevei 



1994 ) or radiation pressure ( Iping & Petterson 1990 : IPringle 



1996). If our results on the rapid realignment of highly twisted 



discs are correct, then this implies that a strong tilting mech- 
anism which operates on short timescales will be required to 
maintain a misaligned disc which can be observed to precess. 

Finally, we briefly consider the implications of our results 
for the formation of planets in misaligned binary systems. The 
early stages of planet formation are believed to involve the 
growth of planetesimals via mutual collisions. As the plan- 
etesimals grow in size, their interaction with the disc through 
gas drag forces will decrease. Planetesimals orbiting at differ- 
ent radii will have their own precession frequencies, likely to 
be different to that of the disc. Consequently, the growth of 
the planetesimals will eventually cause them to develop orbits 
which are increasingly inclined to the disc midplane. The effect 
of this is likely to be an increase in relative motion between 
planetesimals, which will affect the outcomes of mutual colli- 
sions, and the rate at which planetary growth proceeds. A de- 
tailed examination of these issues will be presented in a forth- 
coming paper. 
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